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Abstract. Density profiles of simulated galaxy cluster-sized dark matter haloes are analysed in the context of 
a recently introduced nonextensive theory of dark matter and gas density distributions. Nonextensive statistics 
accounts for long-range interactions in gravitationally coupled systems and is derived from the fundamental 
concept of entropy generalisation. The simulated profiles are determined down to radii of « 1% of R2oo. The 
general trend of the relaxed, spherically averaged profiles is accurately reproduced by the theory. For the main 
free parameter k, measuring the degree of coupling within the system, and linked to physical quantities as the 
heat capacity and the polytropic index of the self-gravitating ensembles, we find a value of —15. The significant 
advantage over empirical fitting functions is provided by the physical content of the nonextensive approach. 
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1. Introduction 

The radial density profiles of dark matter (DM) haloes, 
generated in the process of hierarchical structure forma- 
tion, were modelled in the past primarily on phenomeno- 
logical grounds. High-resolution N-body simulations were 
performed (e.g. Navarro et al. 11996 1997, Moore et al. 
I1999|l in order to reveal the DM density distribution. The 
functional dependence puM {r/rs)~°'{l + r/rs)^'-^^"^ 
fZhao [n)96|l provides reasonable fits to haloes on all ob- 
servable scales (e.g. Ricotti I2003f) . Here rg is a scal- 
ing radius and the parameter a is related to the initial 
power spectrum. For a — 1 this expression represents 
the Navarro, Frenk & White (11^ [TWTjl . hereafter NFW 
profile, while the Moore et al. H1999|l profile follows with 
a = 1.5. Both on theoretical and on observational grounds 
there is still profound discussion on the innermost slope 
of the DM density profile. While numerical simulations 
generally predict a divergent, 'cuspy' behaviour towards 
the centre (e.g. Navarro et al. l2004|l observations favour a 
flat, cored profile (e.g. Gentile et al. I2UU4I and references 
therein). 

Only a few attempts to provide theoretically motivated 
models for density profiles of cosmic structure were pro- 
posed so far. Early analytic work was pioneered by Gunn 
& Gott 1)1972(1 and subsequently elaborated by Fillmore 
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& Goldreich H1984II and Bertschinger H1985|l . who derived 
similarity solutions for the self-similar collapse of spherical 
perturbations. Those solutions were naturally power-laws 
for virialised objects. Later, correlations between the slope 
of the density profile and the form of the power spectrum 
of the initial density perturbations were recognized (e.g. 
Hoffman I1988|l . More recent attempts based on analytic 
infall models were provided by Williams et al. fZ004) and 
Ascasibar et al. (|2004|l while Hansen (2004,) and Austin et 
al. H2005|l used the Jeans equation to gain further insight. 

In this work we apply a recently proposed nonexten- 
sive theory for DM and gas density profiles (Leubner 
[2005)) , accounting for long-range interactions and corre- 
lations present in any astrophysical system, to results 
obtained from N-body simulations. The term 'nonexten- 
sive' is used for statistical ensembles, whose total entropy 
is not an additive quantity as in standard Boltzmann- 
Gibbs-Shannon (BGS) thermo-statistics. Instead, pseudo- 
additive entropy generalisation accounts for correlations 
between sub-systems. The possibility to describe DM den- 
sity profiles using nonextensive statistics was recently 
also suggested (Hansen et al. 2003ll in the context of 



Eddingtons formula (e.g. Binney & Tremaine I1994|l . 
whereas an earlier theoretical link between entropy gener- 
alisation and self-gravitating structures was provided by 
Plastino & Plastino itPMIl. 



Send offprint requests to: 
Thomas . Kronberger@uibk . ac . at 



The classical BGS extensive thermo-statistics con- 
stitutes a powerful tool when microscopic interactions 
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and memory are short ranged and the environment is 
an EucHdean space-time, a continuous and differentiable 
manifold. However, in the present situation we are deal- 
ing with astrophysical systems, generally subject to spatial 
or temporal long-range interactions, i.e. ensembles evolv- 
ing under correlations that makes their behavior nonex- 
tensive. A suitable generalization of the BGS entropy for 
statistical equilibrium was first proposed by Renyi (1955 ) 
and subsequently by Tsallis (1988), preserving the usual 
properties of positivity, equiprobability and irreversibility, 
but suitably extending the standard extensivity or addi- 
tivity to nonextensivity. The main theorems of the clas- 
sical Maxwell-Boltzmann statistics admit profound gen- 
eralizations within nonextensive statistics (sometimes re- 
ferred to as q-statistics where q characterizes the degree 
of nonextensivity of the system), wherefore a variety of 
subsequent analyses were devoted to clarify the mathe- 
matical and physical consequences of pseudo-additivity, 
for an early review see e.g. Tsallis (|1995f) . Those include a 
reformulation of the classical N-body problem within the 
extended statistical mechanics (Plastino et al. 1994) and 
the development of nonextensive distributions (Silva et al. 
119981 Almeida I2001|l . A deterministic connection between 
the generalized entropy and the resulting power-law func- 
tionals was recognized (Andrade et al. 2002), as well as the 
dual structure of nonextensive statistical theory (Karlin et 
al. lMi^ . 

Besides new insights into fundamental physics as- 
trophysical applications support q-non-extensive statis- 
tics, proposing distribution functions for stellar polytropes 
(Plastino et al. .1993,) and explaining the solar neutrino 
counting rate (Kaniadakis et al. ll996|) . Cosmological tests 
provide a bound on the degree of possible non-extensivity 
from primordial helium abundance data (Torres et al. 
I1997|l and the temperature of the cosmic microwave back- 
ground radiation in a Robertson- Walker universe was 
shown to be independent of the degree of non-extensivity 
(Hamity et al. I1996f) . Furthermore, generalized statistics 
was applied to large scale astrophysical systems subject 
to long-range gravitational forces in view of galaxy dis- 
tribution studies (Nakamichi et al. I2U(J2|) and reformu- 
lated also in the context of special relativity (Kaniadakis 
12002 ). Recently nonextensive theory was successfully in- 
troduced also to study the scale dependence of intermit- 
tent flows in astrophysical plasma turbulence, appearing 
as consequence of long-range interactions (Leubner and 
Voros I2005|l . In this context the resulting bi-kappa dis- 
tribution function appears as manifestation of the dual 
nature of nonextensive statistics, which provides also the 
physical background of entropy bifurcation in the theoret- 
ical context of DM and gas density distributions of clus- 
tered matter fLeubner l2005|) . 

We relate in the following nonlocality in gravitationally 
clustered astrophysical structures to the presence of long- 
range forces in nonextensive systems and demonstrate 
that density distributions derived within the framework 
of entropy generalisation consistently reproduce the den- 
sity profile of cluster-sized dark matter haloes forming in 



A cold dark matter (ACDM) N-body simulations. This re- 
sult is furthermore supported by the integrated mass pro- 
file of an observed relaxed galaxy cluster. Consequently, 
nonextensive statistics provides a physically interpretable 
alternative to empirical fitting procedures. 

2. Theory 

Let us first illuminate the property of pseudo-additivity 
in the context of nonextensive entropy generalization by 
considering two sub-systems A and B such that 

S^{A + B)= S^{A) + S^{B) + -S^{A)S^{B) (1) 

where denotes the entropy as depending on the 
entropic index k, a parameter specifying the degree of 
nonextensivity in the system. For k — oo the last term 
on the right hand side cancels leaving the additive en- 
tropy of the standard BGS statistics. Hence, nonlocality or 
long-range interactions are introduced by the multiplica- 
tive last term on the right hand side of Eq. (QJ account- 
ing for correlations between the subsystems A and B. As 
measure of the entropy mixing the parameter k quanti- 
fies the degree of statistical correlations in the system and 
thus accounts for nonlocality and long-range interactions 
or couplings, respectively. In general, the pseudo- additive, 
K-weighted term may assume positive or negative defi- 
nite values indicating a nonextensive entropy bifurcation. 
Obviousely, nonextensive systems are subject to a dual 
nature since positive K-values imply the tendency to less 
organized states where the entropy increases whereas neg- 
ativ K- values provide a higher organized state of decreased 
entropy, see Leubner ('2005 ). 

Next, consider a DM halo as a self-gravitating coUi- 
sionless system of particles in dynamical equilibrium (e.g. 
Burkert Onnni Firmani et al. IMlUj) . Consistent with Eq. 
(^) the corresponding generalised entropy S{k), charac- 
terising systems subject to long-range interactions and 
couplings in nonextensive statistics, reads (Tsallis 119881 
Leubner 12004)1 

S,^KkB{J2pl'^''~l) (2) 

where pi is the probability of the i*^ microstate, fc^ is 
Boltzmann's constant. The transformation k = 1/(1 — q) 
links the K-formalism, commonly applied in astrophysi- 
cal plasma modelling to the Tsallis q-statistics (Leubner 
'2002') . Here k — oo represents the extensive limit of statis- 
tical independence and recovers the classical BGS entropy 

as Sb = -ksJ^Pi^'^Pi- 

Since entropy and probability distributions reside 
physically on the same level the corresponding general- 
ized energy distribution, accounting for long-range inter- 
actions, is available. In Maxwells derivation the veloc- 
ity components of the distribution /(v) are uncorrelated 
where Inf can be expressed as a sum of the logarithms of 
the one dimensional distribution functions. In nonexten- 
sive systems one needs to keep correlations between the 
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components, which can be done conveniently by extrem- 
izing the entropy under conservation of mass and energy 
yielding the corresponding one- and three-dimensional 
power-law distributions in velocity space fLeubner I2(JU4|I 
Here we must retain the spatial dependence and con- 
sider a spherical symmetric, self-gravitating and collision- 
less N-body system where the corresponding steady state 
phase-space distribution f{r,v) obeys the Vlasov equa- 
tion. If the system of particles itself provides the gravi- 
tational potential $ and /(r, v) is regarded as the mass 
distribution then Poisson's equation reads 



A$ = AttGp = AttG 



(3) 



This representation governs the equilibrium of the sys- 
tem where commonly the relative particle energy Er — 
— l/2w^ -I- ^E" is introduced (e.g. Binney & Tremaine .l994j) . 

Extremizing the generalised entropy functional ^ 
with regard to conservation of mass and energy (Plastino 
the resulting energy distribution f^{Er) involved in 
Eq. © reads 



1 



1 - * 



(4) 



The superscripts refer to the positive or negative in- 
tervals of the entropic index k, accounting for less 
and higher (-) organized states and thus reflecting the 
accompanying entropy increase or decrease, respectively, 
(see Leubner I2005II . a represents the mean energy or vari- 
ance of the distribution and denote the corresponding 
normalisation constants depending on the entropic index 
K, for details see Leubner (21)1) 1"). For k ^ oo Eq. Q 
approaches the exponential distribution function defining 
the density profile of the isothermal sphere (Binney & 



Tremaine ll994|l . 

After incorporating the sign of k into Eq. Q , we per- 
form separately for positive and negative definite k the 
integration over all velocities where B^ must be used con- 
sistently. The resulting solution provides the relation for 
the density evolution of a system subject to long range 
interactions in a gravitational potential as 



= Po 



1* 



3/2-K. 



(5) 



Eq. (01 generates for finite positive values of k pro- 
nounced density tails, whereas for negative K-values the 
solutions are restricted within the cutoff at k = 5* /ct^ and 
K = —oo. 

The duality of equilibria in nonextensive statistics is 
manifest in two families, the nonextensive thermodynamic 
equilibria and the equilibria of kinetic equations, where 
positive K-values correspond to the stationary states of 
thermodynamics and negative K-values to kinetic station- 
ary states (Karhn et al. I2002|l . The Hmiting BGS state 
for K = oo is characterized by self-duality. Physically 
the parameter k is related to the heat capacity of the 
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Fig. 1. Characteristics of the radial DM density profile: A 
sequence of solutions to Eq. © is plotted corresponding 
to different values of the parameter k, varied from -2 to 
— oo (left to right). The second free parameter a was kept 
fixed for this schematic plot. 



medium where negative heat capacity, corresponding to 
values K < 0, is a typical property of self-gravitating sys- 
tems, see e.g. Firmani et al. ()2(J(J(Jfl . The nonextensive 
bifurcation into two distributions f^{Er) or p^, respec- 
tively, requires to identify the density profile (O for pos- 
itive definite k as the proper distribution of the thermo- 
dynamic state of the gas, whereas the negative definite 
counterpart is associated with the self-gravitating DM dis- 
tribution, we will focus on. In the limit k — > cx) both solu- 
tions merge at the isothermal sphere solution, defined by 
the limiting exponential distribution obtained for k — oo 
in Eq. gj). 

Finally, we combine Poisson's equation A^f = —AnGp 
and the density distribution Eq. (|3J), which yields after 
re-arranging for ^' a second order nonlinear differential 
equation for the variation of a spherically symmetric gas 
and DM distribution as fLeubner I2()()5|l 



dr'^ r dr 



n p dr 



A-kGu 



(3/2 - n)a' 



Po 



0(6) 



where n — 3/2 — k is introduced and corresponds to the 
polytropic index of stellar dynamical systems. For nega- 
tive definite K-values the DM density profiles defined by 
Eq. @ are found numerically by solving the corresponding 
two first order differential equations with a Runge-Kutta 
routine. Fig.^shows the characteristics of the radial den- 
sity profile as depending on the variation of the entropic 
index n and computed from Eq. jnj for DM profiles. With 
increasing k, the solution converges from left to right to the 
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isothermal sphere density profile, an environment subject 
to statistical independence of the constituents. 

3. N-body simulations 

The numerical simulations were performed using an 
adapted version of the treecode of Barnes & Hut H1986|l , 
using 128'^ particles in a comoving volume of (32 h~^ 
Mpc)^ in a ACDM cosmology with Ha = 0.7, fim = 
0.3, h = 0.7, CT8 = 0.93. We employ the constrained ran- 
dom field technique of Hoffman & Ribak l|1991|l . as im- 
plemented by van de Weygaert & Bertschinger H1996|l . to 
form a rich cluster in the centre of the simulation box. 
The (Plummer) softening length was set to 7 h~^ kpc, 
which also sets the spatial resolution of the simulation. 
Well over a thousand time steps were used to integrate 
the constrained initial conditions to the present epoch. 
The diagonal components of the velocity dispersion tensor 
for the mass within 3 Mpc are: ax = 804 km/s, ay ~ 872 
km/s, and az = 1043 km/s. 

To obtain a spherically averaged mass density profile 
at each time step, we first determine the cluster centre 
using an iterative technique similar to the one presented 
in Power et al. H2003|l . We compute the centre of mass from 
particles within spheres of decreasing radius. Then we bin 
the particles logarithmically according to their distance 
from the cluster centre and calculate the density profile 
out to a radius of ~ 2.5 Mpc (physical distance). In order 
to get reliable results, we consider only scales larger than 
the softening length of the particles and the number of 
particles within each bin is always well above 100. 

4. Discussion and conclusions 

With this setup we studied a sample of 3 cluster-sized DM 
haloes with different formation histories. We confirm the 
commonly known property that the inner density profile of 
the halo does not evolve anymore after a redshift of about 
one. We then use the relaxed density profile of the most 
massive cluster halo of our sample (Mt,ir=8.2x 10^"* M0 
within R200 at z=0), which is also the most virialised sys- 
tem already at zwl, and fit the noncxtcnsive distribution 
by solving the differential equation 10 . The result is shown 
in Fig. 12 together with the corresponding NFW profile for 
comparison. The nonextensive density profile follows the 
result of the simulation accurately at all cluster radii with 
best fitting values of k = —15, indicating the presence of 
strong correlations in the system, and a = 0.12. Hansen 
et al. H2005|l investigated possible theoretical constraints 
on the entropic index q, or equivalently k. Arguing that 
negative velocity dispersions are difficult to interpret in a 
cosmological context, he arrived at the condition q < 5/3, 
which reads in the k formalism k < —3/2. By demanding 
a positive proportion between temperature and internal 
energy Boghosian (1999) arrived at the same constraint. 
Our best fitting value of k = —15 clearly obeys this re- 
striction. 
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Fig. 2. Radial dark matter density profile obtained from 
the simulation (crosses). The solid line shows a fit of the 
nonextensive theory to the data with best fitting values 
of K = —15 and a = 0.12. For comparison, also the best 
fitting NFW profile is provided as a dashed line but shifted 
to the right for better visibility. i?200 indicates the virial 
radius. 



Fig-Eldemonstrates that also the empirical NFW pro- 
file represents the simulated radial density distribution 
well. However, the use of the nonextensive approach al- 
lows one to interpret differences that might occur, e.g. in 
comparison to observations, within the physical context of 
entropy generalisation. On the other hand, in deriving Eq. 
(jSJ we assumed that the phase-space distribution function 
/(r, v) is a function of energy only. This is essentially true 
for systems having an isotropic velocity dispersion tensor. 
If this is not the case, as suggested by numerous numerical 
simulations (e.g. Cole & LaceyEMI Natarajan, Hjorth & 
van Kampen I1997|l . f{r,v) is generally also a function of 
the angular momentum vector L of the system. Therefore, 
for systems with a highly anisotropic velocity dispersion, 
solutions of Eq. JHJ might not represent the density pro- 
file in detail. Furthermore, there is no physical reason to 
keep K constant, as used within the present analysis. For 
instance, one might suspect k to vary with radius, as the 
strength of the correlations could be a function of the dis- 
tance. However, since the shape of the density profiles, 
from simulations and observations, are reproduced accu- 
rately within the applied simplifications (see Figs. El and 
13), and given the fact that also the cluster halo shown in 
Fig.|2]is not completely isotropic (see section 4), modifica- 
tions due to these open issues are expected to be small. In 
addition, presently the resolution of the simulation is in- 
sufficient to obtain robust estimates for the density distri- 
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Fig. 3. Integrated mass profile of the galaxy cluster A1413 
(crosses, Pointecouteau et al. I2005|l compared to the sim- 
ulated (dashed) and theoretical (solid) curve. 



bution at smaller radii than presented in Fig. |21 However, 
from high resolution simulations (e.g. Navarro et al. l2(JU4|l 
it is known that DM profiles do not become flat at cur- 
rently probed radii (<1% of R2oo)- Finally, in Fig. |31the 
present analysis is tested on an integrated mass profile 
evaluated from X-ray observations. Along with the inte- 
grated simulation profile and the nonextensive solution 
a normalized profile of the galaxy cluster A1413 is plot- 
ted (Pratt & Arnaud 2002; Pointecouteau et al.^i)5J. In 
the observable radial range both, the simulated and the 
nonextensive profile, perfectly match the observed data 
set. Note that A1413 is a relaxed cluster with no cooling 
core, consistent with the theoretical context restricted to 
equilibrium states. 

In principle, as the topology of the solutions to equa- 
tion JHJ is very complex, also density profiles from sys- 
tems in earlier evolutionary stages can be well reproduced. 
However, in deriving the differential equation © we im- 
plied that the system is relaxed, as we extremized the 
generalised entropy functional (j^J to get the equilibrium 
phase space distribution function. As the cluster is not re- 
laxed at redshifts z > 1 the physical interpretation of the 
fitting parameters fails or is extremely complicated in this 
redshift range. 

We also note that the strength of correlations within a 
system depends on the cosmological model, or more pre- 
cisely on the form of the power spectrum of the initial 
density fluctuations, i.e. P(fc) oc fc". It has been recog- 
nized already in early numerical works that increasing n 
leads to a steepening of the density profiles (e.g. Hoffman 



119881 Efstathiou et al. I1988|l . Our theoretical framework 
of non-extensive statistics and the derivation of Eq. lO 
do not depend on the initial power spectrum. Therefore 
the theory is also valid for different power spectra than 
the one adopted for this work {n= 1). Even without per- 
forming additional simulations ourselves, we can conclude 
that long-range interactions are more important for high 
spectral indices n, as steeper profiles are realised by lower 
values of k. 

In conclusion, nonextensive statistics provides on phys- 
ical grounds access to the study of DM density profiles in 
relaxed clusters and is able to model the corresponding 
equilibrium states of self gravitating coUisionless systems. 
The significant advantage over empirical fitting functions 
is provided through the physical content of the parameters 
involved in the nonextensive approach, with k as a mea- 
sure of the degree of correlations in the system and cr, the 
characteristic energy. Links between the entropic index k, 
the heat capacity and the polytropic index of DM haloes 
are subject of a further study. 
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